//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : solv_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef SOLV_KERNELS_H
#define	SOLV_KERNELS_H

namespace rrbsolver
{

/* ===========================================================================
// In this routine the r2 nodes are updated via:
//
//                        r1+ld      r1+ld+1
//                             \    /
//        (rather than se) g.ee \  / g.ss (rather than sw)
//                               r2
//        (rather than ne) g.nn /  \ g.ww (rather than nw)
//                             /    \
//                           r1      r1+1
// 
// For example, r2 is affected by its north-east neigbour via a factor g.ss
// rather than g.sw, because previously (in the preconditioner) we have set
// nn -> ne, ee -> se, ww -> nw, ss -> sw.
//
// The result is written in the next coarser grid (level down). 
// The computation is:
//
// vecZnew[r2] = vecR[r2] - g.ss[r1 + ld + 1] * vecR[r1 + ld + 1]
//                        - g.ww[r1 + 1]      * vecR[r1 + 1]
//                        - g.ee[r1 + ld]     * vecR[r1 + ld]
//                        - g.nn[r1]          * vecR[r1]
//
// In the next code:   x  = 2x2 grid repr. of vector R (incoming residual)
//                     z  = current (fine) 2x2 grid repr. of vector Z
//                     y  = next (coarser) 2x2 grid repr. of vector Z
//                     g  = current (fine) preconditioner grid
//                     gg = next (coarser) preconditioner grid (note: stencils 
//                          not used, but the grid comes in handy for knowing
//                          the next leading dimension (ldd), etc.                            
// =========================================================================*/
template <class T>
__global__ void kernel_solv1v1(T *y, T *z, const T *x, 
    const Grid gg, const Grid g)
{
    int ld = g.nx;
    int ldd = gg.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;

    const T x_r1up1vp1 = __FETCH_XX(u_r1+1, v_r1+1, ld);
    const T x_r1up1    = __FETCH_XX(u_r1+1, v_r1  , ld);
    const T x_r1vp1    = __FETCH_XX(u_r1  , v_r1+1, ld);
    const T x_r1       = __FETCH_XX(u_r1  , v_r1  , ld);

    int loc_new = ldd * (BORDER_WIDTH + by * (By >> 1)) 
                   + (ty & 1) * (BORDER_WIDTH + gg.cy) * ldd 
                   +  BORDER_WIDTH + bx * (Bx >> 1)  
                   + (tx & 1) * (BORDER_WIDTH + gg.cx) 
                   + (ty >> 1) * ldd + (tx >> 1);

    T sum = x[loc_r2] - g.ss[loc_r1 + ld + 1] * x_r1up1vp1  // north-east
                      - g.ww[loc_r1 + 1     ] * x_r1up1     // south-east
                      - g.nn[loc_r1         ] * x_r1        // south-west
                      - g.ee[loc_r1 + ld    ] * x_r1vp1;    // north-west

    // this is part of dividing by diagonal
    z[loc_r1] = x_r1 * g.cc[loc_r1]; 

    y[loc_new] = sum;
}


// ===========================================================================
// This routine is exactly the same as solv::kernel_solv1v1, except from 
// that now is read from the previous 2x2 grid repr. of vector Z
// instead of vector R (like is done in solv::kernel_solv1v1)
// ===========================================================================
template <class T>
__global__ void kernel_solv1v2(T *y, const T *x, const Grid gg, const Grid g)
{
    int ld = g.nx;
    int ldd = gg.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;

    const T x_r1up1vp1 = __FETCH_XX(u_r1+1, v_r1+1, ld);
    const T x_r1up1    = __FETCH_XX(u_r1+1, v_r1  , ld);
    const T x_r1vp1    = __FETCH_XX(u_r1  , v_r1+1, ld);
    const T x_r1       = __FETCH_XX(u_r1  , v_r1  , ld);

    int loc_new = ldd * (BORDER_WIDTH + by * (By >> 1)) 
                   + (ty & 1) * (BORDER_WIDTH + gg.cy) * ldd 
                   +  BORDER_WIDTH + bx * (Bx >> 1)  
                   + (tx & 1) * (BORDER_WIDTH + gg.cx) 
                   + (ty >> 1) * ldd + (tx >> 1);

    T sum = x[loc_r2] - g.ss[loc_r1 + ld + 1] * x_r1up1vp1  // north-east
                      - g.ww[loc_r1 + 1     ] * x_r1up1     // south-east
                      - g.nn[loc_r1         ] * x_r1        // south-west
                      - g.ee[loc_r1 + ld    ] * x_r1vp1;    // north-west

    y[loc_new] = sum;
}


// ===========================================================================
// This routine is exactly the same as solv::solv_ph1v1, except from that now
// no next grid is filled. This routine is thus used right before the final
// level solver is called.
// ===========================================================================
template <class T>
__global__ void kernel_solv1v3(T *x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;

    const T x_r1up1vp1 = __FETCH_XX(u_r1+1, v_r1+1, ld);
    const T x_r1up1    = __FETCH_XX(u_r1+1, v_r1  , ld);
    const T x_r1vp1    = __FETCH_XX(u_r1  , v_r1+1, ld);
    const T x_r1       = __FETCH_XX(u_r1  , v_r1  , ld);

    T sum = x[loc_r2] - g.ss[loc_r1 + ld + 1] * x_r1up1vp1  // north-east
                      - g.ww[loc_r1 + 1     ] * x_r1up1     // south-east
                      - g.nn[loc_r1         ] * x_r1        // south-west
                      - g.ee[loc_r1 + ld    ] * x_r1vp1;    // north-west

    x[loc_r2] = sum;
}


// ===========================================================================
// In this routine the red nodes are updated via:
//
//                   b2                              b1+ld
//                   |                                 |                     
//                   | g.nn                            | g.nn
//                   |                                 |  
//       b1-1 ______ r1 ______ b1     and    b2 ______ r2 ______ b2+1
//            g.ww   |    g.ee                  g.ww   |    g.ee
//                   |                                 |
//                   | g.ss                            | g.ss
//                 b2-ld                               b1
//
// or, together,
//
//                                      b1+ld
//                                        |
//                                        | g.nn
//                                        |
//                              b2 ______ r2 ______ b2+1
//                              |  g.ww   |    g.ee
//                         g.nn |         |
//                              |         | g.ss
//                  b1-1 ______ r1 ______ b1
//                       g.ww   |    g.ee
//                              |
//                              | g.ss
//                            b2-ld
//
// For example, r2 is affected by its north neighbour via a factor g.nn.
// Note that the left neigbhour of r2 is the north neighbour of r1, namely b2,
// likewise, the east neigbbour of r1 is the south neighbour of r2, namely b1.
//
// The result is written in the same grid as is read from. The computation is:
//
// for r1 nodes:
//
// vecZ[r1] = vecZ[r1] - g.ss[b2]      * vecZ[b2]
//                     - g.ww[b1]      * vecZ[b1]
//                     - g.nn[b1 - ld] * vecZ[b2 - ld]
//                     - g.ee[b2 - 1]  * vecZ[b2 - 1]
// 
// for r2 nodes:
//
// vecZ[r2] = vecZ[r2] - g.ss[b1 + ld] * vecZ[b1 + ld]
//                     - g.ww[b2 + 1]  * vecZ[b2 + 1]
//                     - g.nn[b1]      * vecZ[b1]
//                     - g.ee[b2]      * vecZ[b2]
// 
// In the next code:   x  = 2x2 grid repr. of vector Z
//                     g  = preconditioner grid
// ===========================================================================
template <class T>
__global__ void kernel_solv2(T *x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;
    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;
    int loc_b1 = ld * v_b1 + u_b1;
    int loc_b2 = ld * v_b2 + u_b2;

    const T x_b1    = __FETCH_XX(u_b1,   v_b1, ld);
    const T x_b1um1 = __FETCH_XX(u_b1-1, v_b1, ld);
    const T x_b2    = __FETCH_XX(u_b2,   v_b2, ld);
    const T x_b2up1 = __FETCH_XX(u_b2+1, v_b2, ld);
    const T x_b2vm1 = __FETCH_XX(u_b2, v_b2-1, ld);
    const T x_b1vp1 = __FETCH_XX(u_b1, v_b1+1, ld);

    x[loc_r1] -= g.ss[loc_b2]      * x_b2 +    // north
                 g.ww[loc_b1]      * x_b1 +    // east
                 g.nn[loc_b2 - ld] * x_b2vm1 + // south
                 g.ee[loc_b1 - 1 ] * x_b1um1;  // west

    x[loc_r2] -= g.ss[loc_b1 + ld] * x_b1vp1 + // north
                 g.ww[loc_b2 + 1 ] * x_b2up1 + // east
                 g.nn[loc_b1]      * x_b1 +    // south
                 g.ee[loc_b2]      * x_b2;     // west
}


// ===========================================================================
// Note: this routine is the backward counterpart of the forward routine 
// solv::kernel_solv2()
//
// In this routine the black nodes are updated via:
//
//                   r2                               r1+ld
//                   |                                  |                     
//                   | g.nn                             | g.nn
//                   |                                  |  
//         r1 ______ b1 ______ r1+1   and   r2-1 ______ b2 ______ r2
//            g.ww   |    g.ee                   g.ww   |    g.ee
//                   |                                  |
//                   | g.ss                             | g.ss
//                 r2-ld                                r1
//
// or, together,
//
//                                r1+ld
//                                  |
//                                  | g.nn
//                                  |
//                      r2-1 ______ b2 ______ r2
//                           g.ww   |    g.ee |
//                                  |         | g.nn
//                                  | g.ss    |
//                                  r1 ______ b1 ______ r1+1
//                                     g.ww   |    g.ee
//                                            |
//                                            | g.ss
//                                          r2-ld 
//
// For example, b2 is affected by its north neighbour via a factor g.nn.
// Note that the right neigbhour of b2 is the north neighbour of b1, namely r2,
// likewise, the west neigbbour of b1 is the south neighbour of b2, namely r1.
//
// The result is written in the next finer grid (level up). The computation is:
//
// for b1 nodes:
//
// vecZ[b1] = vecZ[b1] - g.nn[b1] * vecZ[r2]
//                     - g.ee[b1] * vecZ[r1 + 1]
//                     - g.ss[b1] * vecZ[r2 - ld]
//                     - g.ww[b1] * vecZ[r1]
// 
// for b2 nodes:
//
// vecZ[b2] = vecZ[b2] - g.ss[b2] * vecZ[r1 + ld]
//                     - g.ww[b2] * vecZ[r2]
//                     - g.nn[b2] * vecZ[b1]
//                     - g.ee[b2] * vecZ[r2 - 1]
// 
// In the next code:   x  = 2x2 grid repr. of vector Z
//                     g  = preconditioner grid
// ===========================================================================
template <class T>
__global__ void kernel_solv3(T *y, const T *x, const Grid gg, const Grid g)
{
    int ld = g.nx;
    int ldd = gg.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;
    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_b1 = ld * v_b1 + u_b1;
    int loc_b2 = ld * v_b2 + u_b2;

    int loc = ldd * (BORDER_WIDTH2 + gg.cy + 2 * by * By) 
                   + BORDER_WIDTH2 + gg.cx + 2 * bx * Bx;

    int loc_new;

    const T x_r1    = __FETCH_XX(u_r1,   v_r1, ld);
    const T x_r2    = __FETCH_XX(u_r2,   v_r2, ld);
    const T x_r1up1 = __FETCH_XX(u_r1+1, v_r1, ld);
    const T x_r2vm1 = __FETCH_XX(u_r2, v_r2-1, ld);
    const T x_r1vp1 = __FETCH_XX(u_r1, v_r1+1, ld);
    const T x_r2um1 = __FETCH_XX(u_r2-1, v_r2, ld);

    T sum;

    // b1 nodes:
    sum = x[loc_b1] - g.nn[loc_b1] * x_r2     // north
                    - g.ee[loc_b1] * x_r1up1  // east
                    - g.ss[loc_b1] * x_r2vm1  // south
                    - g.ww[loc_b1] * x_r1;    // west
    loc_new = loc + ldd * (ty << 1) + (tx << 1) + 1;
    y[loc_new] = sum;

    // b2 nodes:
    sum = x[loc_b2] - g.nn[loc_b2] * x_r1vp1  // north
                    - g.ee[loc_b2] * x_r2     // east
                    - g.ss[loc_b2] * x_r1     // south
                    - g.ww[loc_b2] * x_r2um1; // west
    loc_new = loc + ldd * ((ty << 1) + 1) + (tx << 1);
    y[loc_new] = sum;

    // r1 nodes:
    loc_new = loc + ldd * (ty << 1) + (tx << 1);
    y[loc_new] = x_r1;

    // r2 nodes:
    loc_new = loc + ldd * ((ty << 1) + 1) + (tx << 1) + 1;
    y[loc_new] = x_r2;
}


/* ===========================================================================
// Note: this routine is the backward counterpart of the forward routine 
// solv::solv_ph1()
//
// In this routine the R2 nodes are updated via:
//
//                         r2-1      r2
//                             \    /
//        (rather than se) g.ee \  / g.ss (rather than sw)
//                               r1
//        (rather than ne) g.nn /  \ g.ww (rather than nw)
//                             /    \
//                      r2-ld-1      r2-ld
// 
// For example, r1 is affected by its north-east neigbour via a factor g.ss
// rather than g.sw, because previously (in the preconditioner) we have set
// nn -> ne, ee -> se, ww -> nw, ss -> sw.
//
// Moreover, g.ee[r2-1] (actually south-east) is the same as g.ww[r1] (actually 
// north-west). Therefore, we can use the ee, ss, nn, ww contributions given
// in terms of r1 rather than r2, and hence more reads from global memory will
// be coalesced!
//
// The result is written in the same grid as is read from. The computation is:
//
// vecZnew[r1] = vecZ[r1] - g.ss[r2]          * vecZ[r2]
//                        - g.ww[r2 - ld]     * vecZ[r2 - ld]
//                        - g.ee[r2 - 1]      * vecZ[r2 - 1]
//                        - g.nn[r2 - ld - 1] * vecZ[r2 - ld - 1]
// 
// or, faster, using more coalesced memory reads,
//
// vecZnew[r1] = vecZ[r1] - g.nn[r1] * vecZ[r2]
//                        - g.ee[r1] * vecZ[r2 - ld]
//                        - g.ww[r1] * vecZ[r2 - 1]
//                        - g.ss[r1] * vecZ[r2 - ld - 1]
//
// In the next code:   x  = 2x2 grid repr. of vector Z
//                     g  = preconditioner grid                           
// =========================================================================*/
template <class T>
__global__ void kernel_solv4(T *x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;
    int u_r2 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_r2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;

    const T x_r2       = __FETCH_XX(u_r2,   v_r2,   ld);
    const T x_r2vm1    = __FETCH_XX(u_r2,   v_r2-1, ld);
    const T x_r2um1    = __FETCH_XX(u_r2-1, v_r2,   ld);
    const T x_r2um1vm1 = __FETCH_XX(u_r2-1, v_r2-1, ld);

    x[loc_r1] -= g.nn[loc_r1] * x_r2 +       // north-east
                 g.ee[loc_r1] * x_r2vm1 +    // south-east
                 g.ss[loc_r1] * x_r2um1vm1 + // south-west
                 g.ww[loc_r1] * x_r2um1;     // north-west
}

// ===========================================================================
// In this routine the b1 and b2 nodes are divided by the center preconditioner 
// value. This is part of solving Dq = q.
//
// In the next code:   x  = 2x2 grid repr. of vector Z
//                     g  = preconditioner grid
// ===========================================================================
template <class T>
__global__ void kernel_solv5(T *x, const Grid g)
{
    int ld = g.nx;

    int u_b1 = BORDER_WIDTH2 + g.cx + bx * Bx + tx;
    int v_b1 = BORDER_WIDTH + by * By + ty;
    int u_b2 = BORDER_WIDTH + bx * Bx + tx;
    int v_b2 = BORDER_WIDTH2 + g.cy + by * By + ty;

    int loc_b1 = ld * v_b1 + u_b1;
    int loc_b2 = ld * v_b2 + u_b2;

    x[loc_b1] *= g.cc[loc_b1];
    x[loc_b2] *= g.cc[loc_b2];
}

// ===========================================================================
// In this routine the r1 nodes are divided by the center preconditioner value
// This is part of solving Dq = q.
//
// In the next code:   x  = 2x2 grid repr. of vector Z
//                     g  = preconditioner grid
// ===========================================================================
template <class T>
__global__ void kernel_solv6(T *x, const Grid g)
{
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * Bx + tx;
    int v_r1 = BORDER_WIDTH + by * By + ty;

    int loc_r1 = ld * v_r1 + u_r1;

    x[loc_r1] *= g.cc[loc_r1];
}

// ===========================================================================
// Here the final set of most coarse grids is solved in one go.
// The computations are performed on the 2x2 grid repr. of vector Z
// We start from r1/r2/b1/b2 grids having size 16x16 (so the previous
// phases were performed till a grid of size 32x32 remained).
// The coarse grids have size: 16x16, 8x8, 4x4, 2x2, 1x1.
//
// To understand how this kernel works, see:
// - solv::solv_ph1() (either of v1, v2, v3);
// - solv::solv_ph2();
// - solv::solv_ph3();
// - solv::solv_ph4(),
// because this kernel combines them all
//
// In the next code:   x  = 2x2 grid repr. of vector Z
//                     g  = preconditioner grid
// ===========================================================================
template <class T>
__global__ void kernel_solvfinal(T *x, const Grid g)
{
    int ld = g.nx;

    int loc;
    int st = 1;
    int tw = 2;
    int u, v;

    // Part I) solving Lx = x
    for (int nt = DIM_COMPUTE_BLOCK / 2; nt >= 1; nt >>= 1) 
    {
        // Phase 2: for all red nodes
        // This part corresponds with solv::kernel_solv2()
        if (tx < nt) {
            u = (st - 1) + tw * tx;
            v = (st - 1) + tw * ty;
        } else {
            u = (tw - 1) + tw * (tx - nt);
            v = (tw - 1) + tw * ty;
        }
        loc = ld * (BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + v) 
                  + BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + u;

        if (ty < nt && tx < nt << 1) {
            x[loc] -= g.ss[loc + st * ld] * x[loc + st * ld] +
                      g.ww[loc + st     ] * x[loc + st     ] +
                      g.nn[loc - st * ld] * x[loc - st * ld] +
                      g.ee[loc - st     ] * x[loc - st     ];
        }
        __syncthreads();

        // phase 1: for r2 nodes: forward substitution of r1 into r2 nodes
        // This part corresponds with solv::kernel_solv1()
        u = (tw - 1) + tw * tx;
        v = (tw - 1) + tw * ty;
        loc = ld * (BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + v) 
                  + BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + u;    

        if (ty < nt && tx < nt) {
            x[loc] -= g.ss[loc + st * ld + st] * x[loc + st * ld + st] +
                      g.ww[loc - st * ld + st] * x[loc - st * ld + st] +
                      g.ee[loc + st * ld - st] * x[loc + st * ld - st] +
                      g.nn[loc - st * ld - st] * x[loc - st * ld - st];
        }
        __syncthreads();

        st = tw;
        tw *= 2;
    }


    // Part II) solving Dx = x
    // first half
    loc = ld * (BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + ty) 
              + BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + tx;
    T sum = x[loc] * g.cc[loc];
    x[loc] = sum;
    __syncthreads();

    // second half
    loc += ld * (DIM_COMPUTE_BLOCK / 2);
    sum = x[loc] * g.cc[loc];
    x[loc] = sum;
    __syncthreads();


    // Part III) solving L^Tx = x
    st = DIM_COMPUTE_BLOCK / 2;
    tw = DIM_COMPUTE_BLOCK;

    for (int nt = 1; nt <= DIM_COMPUTE_BLOCK / 2; nt <<= 1) 
    {

        // for r1 nodes: backward substitution of r2 into r1 nodes
        if (tx < nt) {
            u = (st - 1) + tw * tx;
            v = (st - 1) + tw * ty;
        } else {
            u = (tw - 1) + tw * (tx - nt);
            v = (tw - 1) + tw * ty;
        }
        loc = ld * (BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + v) 
                  + BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + u;

        if (ty < nt && tx < nt) {
            x[loc] -= g.nn[loc] * x[loc + st * ld + st] +
                      g.ee[loc] * x[loc - st * ld + st] +
                      g.ss[loc] * x[loc - st * ld - st] +
                      g.ww[loc] * x[loc + st * ld - st];
        }
        __syncthreads();

        // for all black nodes
        if (tx < nt) {
            u = (st - 1) + tw * tx;
            v = (tw - 1) + tw * ty;
        } else {
            u = (tw - 1) + tw * (tx - nt);
            v = (st - 1) + tw * ty;
        }
        loc = ld * (BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + v) 
                  + BORDER_WIDTH2 + DIM_COMPUTE_BLOCK + u;

        if (ty < nt && tx < nt << 1) {
            x[loc] -= g.nn[loc] * x[loc + st * ld] +
                      g.ee[loc] * x[loc + st     ] +
                      g.ss[loc] * x[loc - st * ld] +
                      g.ww[loc] * x[loc - st     ];
        }
        __syncthreads();

        tw = st;
        st /= 2;
    }
}

} // namespace rrbsolver

#endif // SOLV_KERNELS_H

